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Abstract 

A path tracking algorithm that adaptively adjusts precision is pre- 
sented. By adjusting the level of precision in accordance with the numer- 
ical conditioning of the path, the algorithm achieves high reliability with 
less computational cost than would be incurred by raising precision across 
the board. We develop simple rules for adjusting precision and show how 
to integrate these into an algorithm that also adaptively adjusts the step 
size. The behavior of the method is illustrated on several examples arising 
as homotopies for solving systems of polynomial equations. 
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1 Introduction 

Path tracking is the task of tracing out a 1 real dimensional solution curve 
described implicitly by a system of equations, typically n equations in n + 1 
variables, given an initial point on, or close to, the path. This can arise in 
many ways, but our motivation is the solution of systems of polynomials via 
homotopy continuation (see El |SJ [7\ |51 E01 EH ) • I n this method, to find the 
isolated solutions of the system f(z) = for given polynomials / : C™ — > C", 
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one constructs a homotopy, H(z, t), H : C n xC^C" such that H(z, 0) = f(z) 
is the target system to be solved while H(z, 1) is a starting system whose isolated 
solutions are known. There is a well-developed theory on how to construct such 
homotopies to guarantee, with probability one, that every isolated solution of 
f(z) = is the endpoint in the limit as t — > of at least one smooth path Zi(t), 
where H(zi(t), t) = on t e (0, 1] and where 2»(1), i = 1,2,..., are the known 
isolated solutions of H(z, 1) = 0. Similar constructions arise in other contexts, 
where the existence of a path leading to the desired solutions may or may not 
be guaranteed. Even when there is no guarantee, experience shows that in some 
application domains continuation techniques yield solutions more reliably than 
Newton's method, especially when good initial guesses are not available. While 
in our applications the path is just a means of arriving at the endpoint, in other 
applications one may desire to accurately trace out the path itself, such as when 
plotting the response of a mathematical model as one of its parameters is varied. 

The most common path tracking algorithms are predictor-corrector meth- 
ods: from an approximate solution point on the path, a predictor gives a new 
approximate point a given step size along the path, then a corrector brings this 
new point closer to the path. For example, one may use an Euler predictor, 
which steps ahead along the tangent to the path, or a higher order predictor 
that uses several recent points and the derivatives of the homotopy function at 
them to extrapolate to the predicted point. Typically, the prediction is then 
used as the initial point for correction by Newton's method. Since the solution 
set is one-dimensional, an extra constraint is introduced to isolate the target of 
the correction. For general homotopies, a useful constraint is to find where the 
solution path intersects the hyperplanc normal to the last computed tangent 
direction. In the more restrictive setting of polynomial systems, the homotopy 
can be designed such that the paths advance monotonically with t, that is, there 
are no turning points, in which case it is acceptable (and simpler) to perform 
corrections by holding t fixed. The adaptive precision algorithm we discuss here 
is compatible with any of these prediction and correction schemes. 

For good results, the predictor step size must be chosen appropriately. Too 
large a step size may result in a prediction outside the zone of convergence of 
the corrector, while too small a step size means progress is slow and costly. 
Consequently, it has long been recognized that adaptive control of the step size 
is crucial for obtaining good reliability without undue computational cost. 

While step size control is well established, less attention has been paid to 
efficient handling of precision. With wider availability of software packages for 
higher precision arithmetic, along with faster computers to execute the software, 
it becomes interesting to consider how adjustable precision might be deployed 
to improve the performance of path tracking algorithms. The issue at stake 
is analogous to step size control: without enough precision, path tracking will 
fail, but the use of excessive precision is inefficient. To address this tradeoff, 
this paper proposes an algorithm that dynamically adjusts the number of digits 
used in computations according to the evolution of the numerical conditioning 
of the homotopy function. 
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In our primary application of interest, the solution of polynomial systems, 
there are several factors driving the need for higher precision. It is well known 
that high degree polynomials often lead to ill-conditioned problems. When 
treating polynomial systems in several variables, the total degree of the system, 
being the product of the degrees of the individual equations, quickly becomes 
large even for low degree polynomials, which can also lead to ill-conditioning. 
Thus, one driving force is the desire to solve larger systems of higher total de- 
gree. A second motivation is that our systems often have some, or possibly 
many, singular solutions, and thus, the solution paths leading to these solu- 
tions are necessarily ill-conditioned near the end. While endgame methods 
exist for enhancing the accuracy with which such endpoints can be estimated, 
for singularities of high enough multiplicity more precision is required. Finally, 
although the homotopy constructions guarantee, with probability one, that no 
path passes exactly through a singularity before reaching its endpoint, there is 
always a chance that a near singular condition can be encountered. To obtain 
the highest reliability possible, we need to detect this and allocate sufficient 
digits to successfully track past such obstructions. 

The paper is organized as follows. In section 2, we review the behavior of 
Newton's method in floating point, revealing how its accuracy and convergence 
properties depend on precision. In section 3, we discuss path tracking with 
adaptive step size control and identify how it fails when precision is insufficient. 
This leads, in section 4, to a novel technique for path tracking using adaptive 
precision. This new adaptive precision path tracking algorithm has been im- 
plemented in a software package, Bertini, currently under development by the 
authors. Several examples are presented in section 5 to illustrate the usefulness 
of adaptive precision. Finally, in section 6, a few related ideas that would make 
interesting studies are discussed. 

2 Background: Newton's method 

The core numerical process in the path tracker is the corrector, which in our case 
is Newton's method. A good predictor speeds up the path tracker by allowing a 
large step while still supplying an initial guess within the convergence region of 
the corrector. However, it is the loss of convergence that causes path tracking 
to fail. In exact arithmetic, as long as the path remains nonsingular, there 
must be a region surrounding the path within which Newton's method converges 
quadratically. With a small enough step At in t, we can be assured of advancing 
along the path, although possibly very slowly. This holds even if we use only 
a zero-th order predictor, i.e., if the point from the last value t k is used to 
initialize the corrector for the new value t k +i = t k + At. In contrast, in inexact 
floating point arithmetic, the convergence region can disappear, thus halting 
the path tracker. Short of this, an unacceptably slow linear rate of convergence 
might dominate, causing the step size to plummet. It can also happen that the 
corrector converges but to an answer that is outside the desired tolerance. 
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Due to these considerations, an analysis of Newton's method in floating point 
is of interest and will help us derive rules for setting the precision used in our 
path tracker. The following analysis resembles that of |T2]- Let F(z) : C" — > 
C™ be continuously differentiable, and denote its Jacobian matrix of partial 
derivatives as J(z). To solve F(z) = by Newton's method given an initial 
guess zq, one iteratively updates Zi, i = 1,2, ... , as 

Solve J{z l )Az l = -F(zi) for Az i} , . 

Zi+i = Zi + Azi. 

Suppose that we work in floating point with unit roundoff u. In other words, 
if we compute with a mantissa of b bits in binary or with P digits in decimal, 
u = 2~ b = 1CP P . We may consider evaluating the residuals F(zi) in higher 
precision, say u < u. Let F(z) be the floating point output of the procedure 
that evaluates F(z). We assume that there exists a function tp depending on z, 
u, and u such that the error e(z) — F(z) — F(z) obeys 

\\e(z)\\<u\\F(z)\\+i;(z,u,u). (2) 

By definition, at a solution point z», we have F{z sr ) = 0, so it is clear that 
the function tp drives the final error. To determine ip, one must examine the 
function F and the program that implements it. We will give a rough rule of 
thumb later for the systems we treat. 

In solving Eq.^for the correction Azj, there is error in evaluating J(Zi) and 
in solving the linear system. Both errors can be absorbed into an error matrix 
Ei such that the computed correction is 

A Zl = {J(zi) + Eiy^Fizi) + e( Zi )). (3) 

We assume this error is bounded by 

\\Ei\\<S(u\\J{zi)\\+<l>(zi,u)), (4) 

for some constant £ > 1 and positive function </>. We expect the first term 
because of roundoff of the Jacobian, whereas <f> accounts for errors in evaluating 
J that do not vanish even when J does. The constant £ accounts for the 
subsequent growth in the error during the linear solve. 

For simplicity of notation, let v = Zi be the current guess, v = z i+1 the new 
guess after a single iteration, and let be the solution point near v. Also, let's 
use the shorthand notations F — F(v), J — J(v), J* = J(i>*), A — \\v — u*| 
and A = \\v — v*\\. In the next paragraph, we will establish a bound on A in 
terms of A. Whenever A < A, the Newton step successfully reduces the error 
in the estimate of the root v* . 

Since F(v*) = 0, the Taylor series of F(z) at v„ gives 
F(z) = J* ■ (z - w») + H.O.T. 
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where the higher order terms, H.O.T., are quadratic or higher inz-t),. Simi- 
larly, 

J(z) = J* + H.O.T. 

where the higher order terms are linear in z — v*. Consequently, in a ball 
B = {z : \\z — V*\\ < R} centered on v* with v £ B, there exist positive 
constants a and (3 such that 

\\F(z)\\ < || J* || || 2 -u* || +a||z-^|| 2 , (5) 
\\F(z) - J*(z - u,)|| < a\\z - v4 2 (6) 
||J*|| < \\J\\ +0\\z-v,\\, \\J- J.|| <0\\z-v.\\. (7) 

In Newton's method, we solve 

(J + E)d = -(F(v) + e) (8) 

for d and take the step 

v = v + d + e, (9) 

where s is the error in forming the sum. The standard model of round-off error 
in floating point addition 18 gives 

\\e\\<u(\\v\\ + \\d\\)<u(A+\\v4 + \\d\\), (10) 

so subtracting v* from both sides of Eq. [51 we have 

A < ||v - + d\\ + u{A + |H| + (11) 

If J is nonsingular and || J _1 ||||i?|| < 1, then (J + E) is nonsingular, the Newton 
step is well defined, and 

\\{J + E)^\\<K\\J-\ jr = __L__ (12) 

Accordingly, from Eq. [SJwe have 

Nl< ^"1(11^1 + INI). (13) 

Also, after adding (J + E)(v — v*) to both sides of Eq. [S]and simplifying using 
Eqs. EH3 w e have 

||t;-t;.+d||<A'||J- 1 ||(||£||A+(a + i 9)A 2 + || e ||). (14) 

Substituting from Eqs. 1131141 into Eq. ^] and using Eqs. 12141 we obtain the 
bound 

A<K\\J- 1 \\{l + u) 2 {a + P)A 2 + 

(K\\ J _1 ||[(2 + £ + u) (u\\J\\ + 4>)} +u)A + K\\ J-^Kl + u)^ + «||«*||. 

(15) 
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This relation holds as long as \\J < 1, so that the linear solve in the 

Newton step is well-defined, and v is in the ball B, so that Eqs. hold. 

If A < A, the Newton step reduces the error in the approximation of the 
root. In exact arithmetic, we have u — <j) — tp — and K = 1, so A < 
|| J~ 1 \\(a + (3) A 2 . The error contracts if the initial guess is accurate enough so 
that A < 1/(|| J^Wia + /?)). If we also have A < 1/(|| J'^a + £)), it is clear 
that all subsequent iterates are nonsingular and contractive, from which one 
has the well-known result that Newton's method converges quadratically to a 
nonsingular solution for a sufficiently accurate initial guess. One sees that the 
more singular the Jacobian is at the root, the slower the convergence and the 
smaller the convergence zone. 

In floating point arithmetic, we cannot expect the error to converge to zero. 
From Eq. one may expect the error to contract until 

Aw (1 +u)K\\J- 1 \\tJj + u\\v4 w J PC||J- 1 ||V' + u||u*||. (16) 

The second term is the error inherent in representing u* in floating point. The 
first term depends on the accuracy, ipj with which the function is evaluated. 
This can be reduced by using higher precision, u, in the function evaluation per 
Eq. [3 The precision of the Jacobian and the linear solve do not affect the final 
error. 

On the other hand, the precision of the Jacobian and the linear solve do affect 
convergence. Without enough precision, ||,/ _1 ||||.E|| may approach or surpass 1, 
which means that the linear solve may fail due to singularity or may yield such an 
inaccurate step that the error diverges. Notice that || J _1 || ||_E|| = u\\ J _1 || || J|| + 
|| J -1 !!^ = uk + ||J _1 ||0, where k — cond(J). The first term, uk, reflects the 
well-known result that in solving a linear system, floating point roundoff is 
magnified by the condition number of the matrix. 

3 Step length control and failure 

To produce an improved path tracking algorithm, it is useful to first examine 
a standard predictor/corrector algorithm to see why adaptive step length con- 
trol generally succeeds when conditioning is mild and why it may fail when 
conditioning is adverse. 

A simple and effective approach for step length control is to adjust the step 
length up or down according to the success or failure of a complete predic- 
tion/correction cycle. Suppose the homotopy function H(z,t) = defines a 
one-dimensional nonsingular path z(t). We are given a start point approxi- 
mately on the path, zq w z(t<j), an ending value of t, and a tolerance to which 
points on the path are to be found. Then, in brief, a predictor/corrector path 
tracker with adaptive step length control may be constructed as follows. 

Initialize Select: an initial step size, s; the number of corrector iterations 
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allowed per step, N > 1; the step adjustment factor, a € (0, 1); the step 
expansion integer, M > 1; and a minimum step size s m i n . 

Predict Estimate a new point near the path whose distance from the current 
point is the step size. 

Correct Iteratively improve the new path point, constraining its distance from 
the prior path point. Allow at most N iterations to reach the specified 
tolerance. 

On success If the tolerance is achieved: 

• Update the current path point to be the newly found path point. 

• If we have reached the final value of t, exit with success. 

• If there have been M successes in a row, expand the step size by 
s = s/a. 

On failure If the tolerance is not achieved: 

• Cut the step length by s = as. 

• If s < s m i n , exit with failure. 

Loop Go back to Predict. 

The key is to allow only a small number of iterations in the corrector, typ- 
ically only N = 2 or TV = 3. This forces the prediction to stay within a good 
convergence region surrounding the path. If a large number of iterations is al- 
lowed, a bad prediction might ultimately converge, but it may wander first and 
become attracted to a completely different path in the homotopy. Keeping N 
small, the step size adaptation slows down to negotiate sharp turns in the path 
and accelerates whenever the path is relatively straight. Properly implemented, 
this results in a robust and efficient path tracking algorithm. 

We can be a bit more precise. Let us do so by specifically considering 
an Euler predictor with a Newton corrector. Both of these derive from the 
linearized local model of the path. The Taylor series at (zi,t\) is 

OH 3H 
H{z 1 + Az,t 1 + At) = H(z 1 ,t 1 ) + —(z u h)Az + —(z 1 ,t 1 )At + H.O.T. (17) 

oz at 

where the higher order terms, H.O.T., are quadratic or higher in (Az,At). 
Ignoring the higher order terms and setting H(z\ + Az, t\ + At) = 0, one has 
the basic Euler predictor and Newton corrector relations. These are a system of 
n equations in n+ 1 unknowns; as long as the combined matrix [dH/dz dH/dt] 
is rank n, there is a well-defined tangent direction and tracking may proceed. 
The predictor adds a constraint on the length of the step along the tangent, 
whereas corrector steps are constrained to move transverse to the tangent. The 
extra constraints are particularly simple in the case where dH/dz is rank n, for 
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then the path progresses monotonically in t, and the step can be controlled via 
the advance of t. Accordingly, one has a linear system to be solved for Az: 



For prediction, we set At — s, the current step size, and for correction, we set 
At = 0. 

Since the neglected terms are quadratic, the prediction error is order 0(s 2 ). 
Thus, in the case of a failed step, cutting the step size from s to as reduces the 
prediction error by a factor of a 2 . In this way, cuts in the step size quickly reduce 
the prediction error until it is within the convergence region of the corrector. 
With a fcth order predictor, the prediction error scales as a k+1 , potentially 
allowing larger step sizes. In any case, the adaptive approach quickly settles 
to a step size s just small enough so that the corrector converges, while the 
next larger step of s/a fails. With a = 1/2 and M = 5, the step size adapts to 
within a factor of 2 of its optimum, with an approximate overhead of 20% spent 
checking if a larger step size is feasible. 

Failure of path tracking with an adaptive step size can be understood from 
the discussion of Newton's method in § El F° r small enough initial error and 
infinite-precision arithmetic, the Newton corrector gives quadratic convergence 
to a nonsingular root. Near a singularity, || J^T 1 ]] is large, which can lead to a 
small quadratic convergence zone and a slower rate of quadratic convergence. 
Inexact arithmetic can further shrink the convergence zone, degrade the conver- 
gence rate from quadratic to linear, and introduce error into the final answer. 
From these considerations, we see that there are two ways for the adaptive step 
size path tracker to halt prematurely near a singularity. 

1. The predictor is limited to a tiny step size to keep the initial guess within 
the convergence zone of the corrector. If this is too small, we may exceed 
the allotted computation time for the path. 

2. The path may approach a point where the final error of the corrector is 
as large as the requested path tracking tolerance. 

The first mode of failure can occur even with infinite precision, but degra- 
dation of the convergence properties with too low a precision increases the oc- 
currence of this failure. The second mode of failure is entirely a consequence of 
lack of precision. By allocating enough precision, we can eliminate the second 
mode of failure and reduce the occurrence of the first mode. It is important 
to note that in some applications there is flexibility in the definition of the 
homotopy, which can be used to enlarge convergence zones and thereby speed 
up path tracking. For example, re-scaling of the equations and variables can 
sometimes help. However, such maneuvers are beyond the scope of this paper, 
which concentrates only on tracking the path of a given homotopy. 




(18) 
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4 Adaptive Precision 



The use of high precision can largely eliminate both types of path tracking failure 
identified above. However, high precision arithmetic is expensive, so it must be 
employed judiciously. One might be tempted to rachet precision up or down in 
response to step failures as in the adaptive step size algorithm. This presents the 
difficulty that there is just one stimulus, step failure, and two possible responses, 
cut the step size or increase precision. In the following paragraphs, we outline 
two possible algorithms for adapting both step size and precision. 

4.1 Adapting precision via path re-runs 

The simplest approach to adapting precision, shown in Figure 1, is to run the 
entire path in a fixed precision with adaptive re-runs. That is, if the path track- 
ing fails, one re-runs it in successively higher precision until the whole path is 
tracked successfully or until limits in computing resources force termination. 
The advantage of this approach is that adaptation is completely external to 
the core path tracking routine. Thus, this strategy can be applied to any path 
tracker that enables requests for higher precision. For example, in the poly- 
nomial domain, the package PHC |16| offers multiple precision, although the 
precision must be set when calling the program. 

The adaptation algorithm of Figure 1 has two main disadvantages. First, 
when too low a precision is specified, the tracker may waste a lot of computa- 
tion near the point of failure before giving up and initiating a re-run in higher 
precision. Second, the whole path is computed in high precision when it may be 
needed only in a small section of the path, often near the end in the approach 
to a singular solution. A slightly more sophisticated treatment can avoid re- 
computing the segment of the path leading up to the failure point by requesting 
the tracker to return its last successful path point. The re-run in higher precision 
can then be initiated from that point on. 

4.2 Stepwise adaptive precision 

Instead of waiting for the adaptive step size method to fail before initiating 
higher precision, we propose to continuously monitor the conditioning of the 
homotopy to judge the level of precision needed at each step. In this way, the 
computational burden of higher precision is incurred only as needed, adjusting 
up and down as the tracker proceeds, while obtaining superior reliability. 

To decide how much precision is needed, we turn to the analysis of Newton's 
method from § [21 We wish to ensure that the achievable accuracy is within the 
specified tolerance and that convergence is fast enough. 

In what follows, we need to evaluate ||J|| and ||J _1 ||. These do not need 
to be very accurate, as we will always include safety margins in the formulas 
that use them. ||J|| is readily available in the max norm, where we use the 
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Figure 1 : Adapting precision via path re-runs 

maximum magnitude of any of its entries. || J~ 1 \\ is more difficult, as we do 
not wish to compute the full inverse of the matrix. This issue has been widely 
studied in terms of estimating the condition number k = 1 1 J 1 1 1 1 J - 1 1 1 . A relatively 
inexpensive method, suggested in |17| and elsewhere, is to choose a unit vector 
b at random and solve Jy = b for y. Then, we use the estimate || J _1 || « \\y\\. 
Although this underestimates || J _1 ||, tests of matrices up to size 10 x 10 show 
the approximation to be reliably within a factor of 10 of the true value, which 
is easily absorbed into our safety margins. 

One requirement is that || J~ 1 \\ \\E\\ should be small enough to ensure that the 
error-perturbed Jacobian is nonsingular. Minimally, we require || J _1 || < 1, 
but by requiring it to be a bit smaller, say || J _1 || \\E\\ < 10 _cri for some a > 1, 
we force K w 1. This removes the growth of K as one possible source of failure. 
Suppose that the error function 4> in Eq. 0] is of the form <fi = $u. Then, our 
first rule is to require 

||J _1 ||f (||J|| < UP" 1 . (19) 

Using P decimal digits of arithmetic results in precision u = 10~ p , so we may 
restate this rule as 

P>a 1 +log 10 [||J- 1 ||£(||J||+4>)]. (A) 

A second requirement is that the corrector must converge within N itera- 
tions, where we keep N small as in the usual adaptive step size algorithm, typi- 
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cally2or3. Let us say that the tolerance for convergence is A = \\v— v*\\ < 1CP T . 
Recall that in each step of Newton's method, we compute d and take the step 
v = v + d. The best estimate available of the accuracy is A sa ||d||, so we declare 
success when ||d|| < 10~ T . Suppose that after i < N iterations this is not yet 
satisfied. We still have N — i iterations to meet the tolerance, and we would like 
to be sure that a lack of precision does not prevent success. Pessimistically, we 
assume that the linear factor in A in Eq. 1151 dominates the quadratic one and 
that the rate of convergence does not improve with subsequent iterations. We 
force K s» 1, and we have tiCl. Including the same safety margin as before, 
lCF 1 , the requirement becomes 

(lO ffl ||J- 1 ||(2 + 5)( U ||J|| +0) + u) JV ~ l ||d|| < l(T r . (21) 

As before, let's assume <f> = $u. Taking logarithms, the number of decimal 
digits of precision must satisfy 

P > a x + log 10 (|| J- 1 ||(2 + £)(\\J\\ + $) + 1) + (r + log 10 ||d||)/(JV - i). (B) 

Since we only apply this formula when the tolerance is not yet satisfied, we have 
||d|| > 10~ r , or equivalently, r + log 10 > 0. This implies that between cor- 
rector iterations, requirement [E] is always more stringent than Eq.^Q However, 
we still use Eq. outside the corrector, because ||d|| is not then available. 

Our third requirement is that the precision must be high enough to ensure 
that the final accuracy of the corrector is within the tolerance at full convergence. 
For this, Eq.^J|is binding, so including a safety margin of 10 _IT2 and using the 
norm of the current approximate solution, \\v\\, as the best available estimate of 
II?;* ||, we require 

||J- 1 ||7/> + u|M| < lO -1 " - " 2 . (23) 

Suppose the error in evaluating the homotopy function is given by tp — ^lu. If 
the function is evaluated in the same precision as the rest of the calculations, 
i.e., u = u, we have the requirement 

P> ( r 2 +T + log 10 (||J- 1 ||*+||i;||). (C) 

If instead we evaluate the function to higher precision, say u = 10~ p < u = 
10 _p , we have the dual criteria 

P>cr 2 +T + log 10 ||i;||, P' > a 2 + T + log 1Q \\J- 1 \\ +log 10 *. (C) 

The effect of adding the two errors is absorbed into the safety factor o^- 

Conditions A, B, and C (or C) allow one to adjust the precision as necessary 
without waiting for the adaptive step size to fail. If necessary, the precision 
can even be increased between corrector iterations. An algorithm using these 
criteria is described by the flowchart in Figure 2. In this flowchart, "Failure" 
in the predictor or corrector steps means that the linear solve of Eq. ^] has 
aborted early due to singularity. Using the magnitude of the largest entry in J 
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as || J|| , Gaussian elimination with row pivoting may declare such a failure when 
the magnitude of the largest available pivot is smaller than u£ \\ J\\, for then the 
answer is meaningless. This is more efficient than completing the linear solve 
and checking condition A or B, as these are sure to fail. 

The algorithm does not attempt corrections at t = 0. This is because in 
our applications the target system F(z) = H(z,0) often has singular solutions. 
It is safer to sample the incoming path while it is still nonsingular and predict 
to t = based on these samples. In this situation, it helps to employ a more 
sophisticated predictor than Euler's method. For example, endgames that esti- 
mate the winding number of the root and use it to compute a fractional power 
series can be very effective ^2 El • 

4.3 Error estimates 

To use the foregoing procedures, we need the function evaluation error, ip, and 
the errors contributing to E, namely, £ and (p. There is a trade-off between 
using rigorously safe bounds for highest reliability or using less stringent figures 
reflecting typical behavior to avoid the overuse of high precision. Rough figures 
are acceptable as this is just a means of setting the precision. Also, a user of the 
path tracker will not usually wish to expend a lot of effort in developing error 
bounds. 

A rigorous and automated way of establishing error bounds is to use interval 
arithmetic. Following that approach, one may wish go all the way and use inter- 
val techniques to obtain a path tracker with fully rigorous step length control, as 
in 0]. However, this can be expensive, due partially to the cost of interval arith- 
metic but more significantly due to the cost of overconservative error bounds, 
which slow the algorithm's progress by driving the step size smaller than nec- 
essary. Still, when rigorous results are desired, it may be worth the cost. The 
method of 0] does not explicitly include adaptive precision, so something along 
the lines discussed here could be useful in modifying that approach. 

Instead of using interval methods, we may approximate errors by accumulat- 
ing their effects across successive operations. Suppose the program to evaluate 
F(z) has been parsed into a straight line program, that is, a sequence of unary 
and binary operations free of branches or loops. Suppose that at some inter- 
mediate stage of computation, we have computed a value a for a real number 
a, such that a lies between (1 — u)a and (1 + u)a. (For a floating point com- 
plex number, this applies to both the real and imaginary parts.) Let's use the 
shorthand a = (1 ± u)a to mean this entire interval. If a = (1 ± u a )a and 
b = (1 ± Ub)b, the product c = ab is computed in floating point with unit 
round-off u as c w ab(l ±max[u, (u a + Ub)]), where the quadratic round-off term 
u a ui, is neglected. The absolute error c — c thus has the approximate bound 
± max[u, (u a + Ub)] \a\ \b\. Similarly, a + b is computed as a + b ± u a a ± u^b which 
has an absolute error bounded by ±max[w|a + b\, {u a \a\ + u b \b\)]. Using just 
these relations, an error bound for any straight-line polynomial function can be 



12 



/Homotopy, start pointA 
\initial constant settings/ 



Check 
(A,C) 




Increase 
precision 



Failure/ Predictor 
Step 



Success 



Violat,eH/y he gf 

v(A,CV, 

'OK 



Failure/ Corrector 
Step 



Violatec 



Success 



Check \QK 
(B,C) 



Consider 
decreasing 
precision 



Consider 
increasing 
Steplength 



Decrease 
Steplength 



Yes 


Failed 




Path 



Yes / Iterations 
^remaining?/ 




Predict 
to t=0 




j^q / Successive 

predictions agree N 



Return 
endpoint 
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calculated in parallel with the function evaluation itself. Similar relations can be 
developed for any smooth elementary function, such as the basic trigonometric 
functions. Assuming that the inputs to the function, including both the input 
variables z and any internal parameters of the function, are all known either 
exactly or with relative round-off error u, the output of the error analysis is 
such that the error in the computed value F{z) is ||F(z)— F{z)\\ = tp(z) — ^(z)u. 
When the result is rounded off to a possibly lower precision u, the total error 
becomes the form shown in Eq. [21 

It is important to note that the error in the function depends on the error 
in its parameters. For example, consider the simple function f(x) = x 2 — 1/3. 
If this is rounded off to g(x) = x 2 — 0.333 before we use high precision to solve 
g(x) — 0, we will obtain an accurate value of V0.333 but we will never get an 
accurate value of Although this is an obvious observation, it can easily be 

forgotten in passing a homotopy function from some application to the adaptive 
precision path tracking algorithm. If coefficients in the function are frozen at 
fixed precision, the algorithm tracks the solutions of the frozen function, not the 
exact problem that was intended. Whether the difference is significant depends 
on the nature of the application and the sensitivity of the function. 

While i\) and <f> concern the errors in evaluating the function and its Jacobian, 
the factor £ concerns the stability of the linear solve. Round-off errors can 
accumulate through each stage of elimination. When Gaussian elimination with 
partial pivoting is used, the worst-case error bound grows as £ — 2" for solving 
an n x n system 0. However, as indicated in 0, £ rarely exceeds n with the 
average case around n~s or . Setting £ = n 2 should therefore be sufficient for 
almost all cases. 

4.4 Application to polynomial systems 

To avoid program complexity and save computation time, it is preferable not to 
perform a full error analysis of the type just described. In many cases a rough 
analysis is sufficient and easily derived. This is indeed possible for the case of 
most interest to us: polynomial systems. 

Suppose h : C n+1 — > C is a degree D homogeneous polynomial 

n 

h(x) = h(xo,xi, . . . ,x n ) = ^c^q 01 • • ■x' 3 X\ y^rfji = D, for all i e I, 

iei j=o 

(26) 

where X is just an index set for the coefficients q. Since h is homogeneous, 
h(Xx) = X D h(x), so if h(x) = 0, then also h(Xx) — 0. Consequently, the solution 
set of a system of homogeneous polynomials can be said to lie in projective 
space P", the set of lines through the origin in C n+1 . Similarly, the solutions of 
multihomogeneous polynomials lie in a cross product of projective spaces, see 

m 

Any inhomogeneous polynomial g(z\, . . . , z n ) can be easily homogenized to 
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obtain a related function G(xo, x\, . . . , x n ), with 

G(l,xi, ...,x n )= g(xi, . . .,x n ). 

Hence, for any solution z» of g(z) — there is a corresponding solution = 
(l,z*) of G{x) = 0. One advantage of homogenization is that we can re-scale 
any solution of G[x) = to make ||x*|| = 1, which often helps numerical 
conditioning. 

Error bounds for homogeneous polynomials can be estimated easily. If we 
rescale x so that the maximum entry in [xq, . . . , x n ) has magnitude 1, then the 
error in evaluating the degree D homogeneous polynomial h(x) as in Eg. 1261 is 
approximately 

iP{x,u)ttuDj2\ci\, i-e., * = £>2M- ( 27 ) 
Similarly, the derivatives have an approximate error bound of 

<l)(x, u) « uD(D i-e-. * = D ( D ~ l ) ( 28 ) 

At first glance, it may seem that errors can be reduced by simply scaling 
the functions, and thereby scaling their coefficients, by some small factor. But 
|| J^T 1 1 will scale oppositely, so the error predicted by Eq. ^|is unchanged. 

5 Computational results 

This section contains a brief discussion of the implementation details for mul- 
tiprecision arithmetic and for evaluating the rules for adapting precision. Then 
we discuss the results of applying the adaptive precision path tracker to three 
example polynomial systems. 

5.1 Implementation details 

Bertini is a software package for computation in numerical algebraic geome- 
try currently under development by the authors and with some early work by 
Christopher Monico. Bertini is written in the C programming language and 
makes use of straight-line programs for the representation, evaluation, and dif- 
ferentiation of polynomials. All the examples discussed here were run using an 
unreleased version of Bertini on an Opteron 250 processor running Linux. 

The adaptation rules, A, B, and C (or C), leave some choices open to the 
final implementation. For the runs reported here, we chose to evaluate function 
residuals to the same precision as the computation of Newton corrections, so rule 
C applied, not rule C. Also, in rules A and B, we chose to use £ — n 2 , where n 
is the number of variables, which is somewhat conservative for typical cases but 
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IEEE single 


IEEE double 


MPFR 


bits 


23 


52 


64 


96 


128 


256 


decimal digits 


6 


15 


19 


28 


38 


77 



Table 1: Number of digits for mantissas at various levels of precision 



underestimates the worst pathological cases. (See Section fOI for more on this 
issue.) The rules require formulas for evaluating the error bounds <fi(x, u) = &u 
and ij)(x, u) = ^>u. These are problem dependent, so we report our choices for 
each of the example problems below. 

To adaptively change precision, Bertini relies on the open source MPFR 
library for multiprecision support. Bertini has data types and functions for 
regular precision (based on the IEEE "double" standard) and higher precision 
(using MPFR). Although the program would be simpler if MPFR data types 
and functions were used exclusively, the standard double precision types and 
functions in C are more efficient, so Bertini uses these whenever the adaptation 
rules indicate that double precision is sufficient. Additional details regarding 
the use of multiple precision may be found using links from the Bertini website. 
Since the use of adaptive precision variables is highly implementation-specific, 
no other details are described here. 

MPFR requires the addition of precision to the mantissa in packets of 32 bits. 
Since the discussion of the examples below involves both binary and decimal 
digits, Table shows how to convert between the two. 

5.2 Behavior of adaptive precision near a singularity 

Consider the polynomial system of Griewank and Osborne 3 , 

( — f§ z i ~~ 2zi02 

This system is interesting because Newton's method diverges for any initial 
guess near the triple root at the origin. A homotopy of the form 

h{z,t) = tg(z) + (l-t)f(z), 

where g(z) is the following start system, constructed as described in 14 : 

((-0.74924187 + 0.13780686i)z 2 + (+0.18480353 - 0.41277609i)ff 2 )* 
((-0.75689854 - 0.14979830i)zi + (-0.85948442 + 0.60841378i)ifi)* 
9 ~ ((+0.63572306 - 0.6281750U>i + (-0.23366512 - 0.46870314i)iT x )* 
((+0.86102153 + 0.27872286i)zi + (-0.29470257 + 0. 336465 78i)iii), 

((+0.35642681 + 0.94511728i)z 2 + (0.61051543 + 0.76031375«)-ff 2 )* 
g2 = ((-0.84353895 + 0.93981958z)zi + (0.57266034 + 0.80575085i)#i)* 
((-0.13349728 - 0.51170231z)zi + (0.42999170 + 0.98290700z)i3i), 
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with Hi and H 2 the extra variables added to the variable groups {zi} and {22}, 
respectively. In the process of homogenizing, two linear polynomials (one per 
variable group) are added to the system and do not depend on t. In this case, 
those polynomials are: 

(-0.42423834 + 0.84693089i)zi +H X - 0.71988539 + 0.59651665i, 
(+0.30408917 + 0.78336869i)z 2 + H 2 + 0.35005211 - 0.52159537i. 

This homotopy has three paths converging to the origin as t — > 0. These paths 
remain nonsingular for t £ (0, 1], so it is interesting to see how the adaptive 
precision algorithm behaves as t approaches zero. 

Using the prescription of Eqs. 1271281 to bound errors, we take D = 3 and 
Sieil c i| ~ 4, so \& = 12 and $ = 24. (This is actually quite conservative, 
because a more realistic bound would be W = 12||z|| and ||z|| is approaching 
zero.) We set the safety digits <j\ = 02 = 0, as these serve only to force precision 
to increase at a slightly larger value of t. 

We set the desired accuracy to r = 8 for t £ [0.1,1] and increased it to 
t = 12 digits thereafter. To watch the behavior of the algorithm for very small 
values of t, we turned off the usual stopping criterion and instead simply ran 
the path to t = 10~ 30 . That is, in the flowchart of Figure 2, after a successful 
correction step, the algorithm was modified to always loop back to step ahead 
in t until t = 10~ 30 . Since for small t and for \\g\\ w 1, the homotopy path has 
\\z\\ sa |i| 1//3 , all three roots heading to the origin were within 10~ 10 of the origin 
at the end of tracking. With an accurate predictor, such as a fractional power 
series endgame the solution at t = can be computed to an accuracy of 
10~ 12 using only double precision, but the purpose of this example is to show 
that the adaptive precision tracking algorithm succeeds beyond the point where 
double precision would fail. 

The result is shown in Figure 3. We plot the right-hand side of rule C, as 
this rule determines the increases in precision for this problem. The jump in 
this value at t = 0.1 is due to our prescribed increase in r at that point. Since 
the path is heading to the origin, if we used the less conservative error estimate 
of ip(z,u) = 12||z||u, increases in precision would be delayed somewhat, as rule 
C would be shifted down by approximately log 10 ||z|| » (l/3)log 10 t. 

In cases where multiple roots are clustered close together or even approaching 
the same endpoint, as in this example, the path tracking tolerance must be kept 
substantially smaller than the distance between roots to avoid jumping between 
paths. Rather than prescribing r versus t, as we did in this example, it would 
be preferable to automatically adjust r as needed. We postpone this for future 
research, but note that the rules given here for adjusting precision should still 
be effective when r is adjusted automatically. 
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Figure 3: Right-hand side of C, condition number, and decimal precision versus \og(t) 



5.3 Behavior of adaptive precision under tighter toler- 
ances 

To illustrate the effect of tightening the tracking tolerance (i.e., increasing r) on 
adaptive precision path tracking, we will consider a polynomial system coming 
from chemistry. To determine chemical equilibria, one may pass from a set of 
reaction and conservation equations to polynomials, as described in 0^]. One 
such polynomial system, discussed in ^2]> is the following: 



/ = 



Uzf + 6z lZ2 + 5zi - 72z| - 18z 2 - 850z 3 + 0.000000002 
0.5zizf + 0.01ziz 2 + 0.13z| + 0.0422 - 40000 
0.03ziz 3 + 0.04z 3 - 850 



As in the previous example, this system was homogenized, although this time 
with only one variable group, so there is just a single homogenizing coordinate. 
Using a compatible total-degree linear product start system, that is, one whose 
polynomials have degrees 2, 3, and 2, respectively, results in a homotopy with 12 
paths. The solution set of this system consists of eight regular finite solutions, 
two of which have some coordinates of size around 10 5 . These eight finite 
solutions are provided in Tabled 

Due to the poor scaling of the problem, the error bounds were set to <!» = 
120, 000 and $ = 240, 000. As opposed to the previous example, the usual 
stopping criterion described in Sectional was employed. No endgame was used, 
though, since the use of endgames speeds convergence. 
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+ 00*i 


z 2 
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+ 
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1 


61788761616c 
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Table 2: The solutions of the chemical system 



Tracking to a final tolerance of 10~ 8 using fixed regular (IEEE double) pre- 
cision, all eight finite solutions were discovered, including the two having large 
coordinates. In tightening the tolerance to 10 -12 , however, the linear algebra 
involved with path tracking broke down for the paths leading to these two large 
solutions. This resulted in step failures and, ultimately, path failure for those 
two paths. However, by increasing precision to 96 bits, all eight regular solutions 
were again found. 

Using adaptive precision with the number of safety digits set to 0, 1, or 2, 
the six finite solutions of moderate size required only one precision increase (to 
64 bits from 52 bits). This increase occurred at the very end of the path. The 
two large finite solutions each needed 96 bits of precision (as expected) . 

Since the initial run in double precision succeeded on the six moderate paths, 
the adaptive method's increase to 64 bits at the end of these paths is not nec- 
essary. For safety, the adaptive rules are designed to be conservative, so some 
extra computational cost is to be expected. 

5.4 A class of univariate polynomials 

The Chebyshev polynomials have been studied extensively and are known to 
have many interesting properties [2]- There is one Chebyshev polynomial in 
each degree, and scaled such that the leading coefficient is 1, they may be 
defined recursively by 

To (a:) := 2, 
T\{x) := x, and 

Ti(x) := xTi-!(x) - T^ 2 (x)/4, for i > 2. 
The solutions of T n+ \(x) are then given by 



cos 



(2n + l - 2fc)?r 
2n + 2 



for k = 0, 1, 
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Degree 


Estimate 


10 


7 


15 


17 


20 


44 


25 


111 


50 


12300 


100 


1.5 • 10 8 


150 


1.9 • 10 12 


200 


2.3 ■ 10 16 


250 


2.8 • 10 20 


300 


3.4- 10 24 



Table 3: Estimates of X^gz \ f° r various degrees 



Several of these polynomials, with degrees ranging from 10 to 300, were 
solved in Bertini using the step-adaptive precision path tracking method devel- 
oped above. The systems were run with a linear homotopy with a total degree 
start system and without homogenizing. In each case, bounds on $ and '5 were 
set almost exactly. In particular, using the estimates of <E> and 'J described in 
Section l4~4l D was chosen exactly, and J^iei l c «l was se * *° the values given in 
Table |3 In each case, both safety digit settings were set to 4 and r was set to 
10. 

The level of precision used by the step-adaptive precision path tracking al- 
gorithm developed above was degree- and path-dependent for the Chebyshcv 
polynomials, although all paths for a given degree needed approximately the 
same level of precision. In each degree that was considered, the path ending 
nearest 1.0 was one of the paths needing the highest level of precision for that 
degree. (This occurs because the spacing between the roots is smallest near 
±1.) The levels of precision used for that path are displayed in Figure 4. It 
should be noted that for every degree considered, a complete solution set with 
all solutions correct to at least 10 digits was found. 

We note that to solve the high degree Chebyshev polynomials, a small initial 
step size was required to get the path tracker started. With too large an initial 
step, the predicted point was so far from the path that the adaptive precision 
rules increased precision to an unreasonable level without ever exiting the cor- 
rector loop. As diagrammed in Figure 2, the algorithm must exit the corrector 
loop before a decrease in step length can be triggered. Various ad hoc schemes 
could detect and recover from this type of error, but we would prefer a step size 
control method based on an analysis of the predictor. For the moment, we defer 
this for future work. 
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Figure 4: Number of bits needed for Chebyshev polynomials of various degrees 

6 Discussion 

On some problems, endgames can speed convergence to the point that singular 
endpoints can be estimated accurately in double precision. It should be noted 
that is not generally the case. Without enough precision, the "endgame operat- 
ing zone" is empty ^3 El- Likewise, endgames based on deflating the system 
to derive a related nonsingular one |S] may need higher than double precision to 
make a correct decision on the rank of the Jacobian at each stage of deflation |3] . 
Moreover, if some other sort of singularity is encountered during path tracking, 
away from t — 0, endgames will not be useful while adaptive precision will be. 
In the case of tight final tolerances or endpoints of paths having high multiplic- 
ity, endgames will again need assistance from higher (and therefore adaptive) 
precision methods. Conversely, high precision is expensive and floating point 
precision can never be made truly infinite, so to get the most out of whatever 
precision one uses, endgames are indispensable. 

The theory going into this new adaptive precision method revolves around 
Newton's method or corrector methods in general. However, corrector methods 
are only one half of basic path tracking. A careful study of predictor methods 
is certainly warranted. The use of different predictor schemes, e.g., Adams- 
Bashforth rather than Euler, is well worth considering. A careful analysis of 
the predictor might be combined with the convergence criteria of the corrector 
to automatically determine a safe step length in place of the trial-and-error 
step length adaptation method we have used here. This might give an efficient 
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alternative to 0], which presents a rigorous step length control algorithm based 
on interval arithmetic. 

Another open question, discussed briefly in Section [5. 21 is the issue of adap- 
tively changing the tracking tolerance in response to close approaches between 
paths. 
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